---
title: "Untitled"
author: "Jacqueline Park"
date: "2024-04-11"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
setwd("/Users/jacqpark/2nd_yr_project")
rm(list=ls())

library(haven)
library(readr)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(desta)
library(maps) # to extract country names
library(countrycode)
library(tradestatistics)
library(WDI)
library(fixest)
library(DBI)
library(RSQLite)
library(states)
library(broom)

load("~/2nd_yr_project/PTD.RData")

```


# PTD
## Data cleansing

```{r}
gsp_df = read_csv("gsp_combined.csv")
names(gsp_df)[names(gsp_df) == "other"] = "other_util_rate"

gsp_cl = gsp_df %>%
  dplyr::group_by(reporter, partner, year) %>%
  dplyr::summarise(total = sum(total),
            dutiable = sum(dutiable),
            gsp_received = sum(gsp_received),
            other_pref = sum(other_pref),
            all_pref = sum(all_pref),
            mfn_rta = sum(mfn_rta),
            mfn_dutyfree = sum(mfn_dutyfree))

## compute PTD
gsp_cl = gsp_cl %>% mutate(
  ptd = gsp_received/total
)

gsp_cl = gsp_cl[!is.na(gsp_cl$ptd),]

## compute 3-yr moving avg of PTD
gsp_cl = gsp_cl %>% 
  group_by(reporter, partner) %>%
  mutate(ptd_ma = zoo::rollmeanr(ptd, 3, fill = NA))

gsp_cl$ptd_ma = round(gsp_cl$ptd_ma, digits = 5)

```



## Merging with other data

```{r}
gsp_cl$iso1 = countrycode(gsp_cl$partner, origin = "country.name", destination = "iso3n")
gsp_cl = gsp_cl %>% select(reporter, partner, iso1, everything())
gsp_cl$reporter = gsub("\\s+", " ", gsp_cl$reporter)

## EU DESTA

eu_desta = desta[desta$entry_type=="base_treaty" &
                  grepl("EC|EU|European Union|European Community", desta$name, ignore.case = F),]
eu_desta = eu_desta %>% filter(!grepl("Economic", name))

# Remove EU duplicates
eu_desta = eu_desta %>% filter(devstat1+devstat2==1)

eu_desta = eu_desta %>% mutate(
    temp = ifelse(devstat1 == 0, country1, country2),
    temp2 = ifelse(devstat1 == 0, iso1, iso2),
    iso1 = ifelse(devstat1 == 0, iso2, iso1),
    iso2 = ifelse(devstat1== 0, temp2, iso2),
    country1 = ifelse(devstat1 == 0, country2, country1),
    country2 = ifelse(devstat1 == 0, temp, country2)
  ) %>%
  select(-c(temp, temp2))

eu_desta = eu_desta %>% distinct(country1, name, .keep_all = T)

eu_desta = eu_desta %>% filter(country1!="Romania" & name!="EC Romania")
eu_desta$country2 = "European Union" 
names(eu_desta)[names(eu_desta)=="country1"] = "partner"
names(eu_desta)[names(eu_desta)=="country2"] = "reporter"

eu_desta = eu_desta %>%
  select(reporter, partner, iso1, year, number, name, contains("ipr_"))

## US DESTA

us_desta = desta %>% filter(country1 == "United States" | country2 == "United States")

us_desta = us_desta %>% mutate(
  country1 = ifelse(country1!="United States", country1, country2),
  iso1 = ifelse(iso1!=840, iso1, iso2),
  country2 = ifelse(country2!="United States", "United States", country2),
  iso2 = ifelse(iso2!=840, 840, iso2)
)

names(us_desta)[names(us_desta)=="country1"] = "partner"
names(us_desta)[names(us_desta)=="country2"] = "reporter"

us_desta = us_desta %>%
  select(reporter, partner, iso1, year, number, name, contains("ipr_"))

both_desta = rbind(eu_desta, us_desta)
both_desta = both_desta %>% select(-partner)

## Merge w/ gsp_cl
gspdat = gsp_cl %>%
  left_join(both_desta, by = c("reporter", "iso1", "year"))

gspdat = gspdat %>%
  left_join(gdppc, by = c("iso1", "year"))

# Identify columns that contain the string "ipr"
ipr_cols <- grep("ipr_", names(gspdat), value = TRUE)

# Replace NA values with zero in the identified columns
gspdat[ipr_cols] <- lapply(gspdat[ipr_cols], function(x) ifelse(is.na(x), 0, x))


```


## PTD prelim analysis

```{r}
## TRIPS
summary(feols(
  ipr_trips_1994_dummy ~ ptd_ma + log(gdppc_ppp+1)|iso1+year,
  data = gspdat, cluster = "iso1"
))
```


```{r}
## Specific enforcement
summary(feols(
  ipr_specific_enforcement_dummy ~ ptd_ma + log(gdppc_ppp+1)|iso1+year,
  data = gspdat, cluster = "iso1"
))

```


```{r}
## Pharma
summary(feols(
  ipr_pharma ~ ptd_ma + log(gdppc_ppp+1)|iso1+year,
  data = gspdat, cluster = "iso1"
))

```

## Throw in covariates

```{r}
## GDP and GDP growth
gspdat = gspdat %>%
  left_join(gdpcons, by = c("iso1", "year"))
gspdat = gspdat %>%
  left_join(gdpgrowth, by = c("iso1", "year"))

gspdat = gspdat %>% select(-iso3c.x, -iso3c.y)
names(gspdat)[names(gspdat) == "NY.GDP.MKTP.PP.KD"] = "gdp"
names(gspdat)[names(gspdat) == "NY.GDP.MKTP.KD.ZG"] = "gdp_growth_pcnt"

## Polity
gspdat = gspdat %>%
  left_join(polity, by = c("iso1", "year"))
gspdat = gspdat %>%
  group_by(partner) %>%
  fill(polity2, .direction = "updown")
gspdat = gspdat %>% mutate(
  democ = ifelse(polity2 > 0, 1, 0)
)


## RTAs signed
gspdat = gspdat %>%
  left_join(desta_dev, by = c("iso1", "year"))

gspdat$rta_signed = replace(gspdat$rta_signed, is.na(gspdat$rta_signed), 0)


```



## PTD analyses w/ covariates

### TRIPS

```{r}
summary(ptdtrips <- feols(
  ipr_trips_1994_dummy ~ ptd_ma +
    log(gdppc_ppp+1) + log(gdp + 1) + gdp_growth_pcnt +
    democ + rta_signed | iso1+year,
  data = gspdat, cluster = "iso1"
))
```



### Specific enforcement

```{r}
summary(ptdenf <- feols(
  ipr_specific_enforcement_dummy ~ ptd_ma +
    log(gdppc_ppp+1) + log(gdp + 1) + gdp_growth_pcnt +
    democ + rta_signed | iso1+year,
  data = gspdat, cluster = "iso1"
))

```



```{r}
summary(ptdspec <- feols(
  ipr_specific_prov ~ ptd_ma +
    log(gdppc_ppp+1) + log(gdp + 1) + gdp_growth_pcnt +
    democ + rta_signed | iso1+year,
  data = gspdat, cluster = "iso1"
))
```



```{r}
summary(ptdsbst <- feols(
  ipr_scope_substantial_dummy ~ ptd_ma +
    log(gdppc_ppp+1) + log(gdp + 1) + gdp_growth_pcnt +
    democ + rta_signed | iso1+year,
  data = gspdat, cluster = "iso1"
))
```




### Pharmaceuticals

```{r}
summary(feols(
  ipr_pharma ~ ptd_ma +
    log(gdppc_ppp+1) + log(gdp + 1) + gdp_growth_pcnt +
    democ + rta_signed | iso1+year,
  data = gspdat, cluster = "iso1"
))
```



```{r}
summary(feols(
  ipr_pharma ~ ptd_ma +
    log(gdppc_ppp+1) + log(gdp + 1) + gdp_growth_pcnt +
    democ + rta_signed,
  data = gspdat
))
```



## PTD - 5-yr moving avg (robustness checks)
```{r}
gspdat = gspdat %>% 
  group_by(reporter, partner) %>%
  mutate(ptd_ma5 = zoo::rollmeanr(ptd, 5, fill = NA),
         ptd_ma7 = zoo::rollmeanr(ptd, 7, fill = NA)) %>%
  dplyr::ungroup()

gspdat$ptd_ma5 = round(gspdat$ptd_ma5, digits = 5)
gspdat$ptd_ma7 = round(gspdat$ptd_ma7, digits = 5)

summary(ptdtrips5 <- feols(
  ipr_trips_1994_dummy ~ ptd_ma5 +
    log(gdppc_ppp+1) + log(gdp + 1) + gdp_growth_pcnt +
    democ + rta_signed | iso1+year,
  data = gspdat, cluster = "iso1"
))

summary(ptdenf5 <- feols(
  ipr_specific_enforcement_dummy ~ ptd_ma5 +
    log(gdppc_ppp+1) + log(gdp + 1) + gdp_growth_pcnt +
    democ + rta_signed | iso1+year,
  data = gspdat, cluster = "iso1"
))

summary(ptdspec5 <- feols(
  ipr_specific_prov ~ ptd_ma5 +
    log(gdppc_ppp+1) + log(gdp + 1) + gdp_growth_pcnt +
    democ + rta_signed | iso1+year,
  data = gspdat, cluster = "iso1"
))

summary(ptdsbst5 <- feols(
  ipr_scope_substantial_dummy ~ ptd_ma5 +
    log(gdppc_ppp+1) + log(gdp + 1) + gdp_growth_pcnt +
    democ + rta_signed | iso1+year,
  data = gspdat, cluster = "iso1"
))

summary(ptdtrips7 <- feols(
  ipr_trips_1994_dummy ~ ptd_ma7 +
    log(gdppc_ppp+1) + log(gdp + 1) + gdp_growth_pcnt +
    democ + rta_signed | iso1+year,
  data = gspdat, cluster = "iso1"
))

summary(ptdenf7 <- feols(
  ipr_specific_enforcement_dummy ~ ptd_ma7 +
    log(gdppc_ppp+1) + log(gdp + 1) + gdp_growth_pcnt +
    democ + rta_signed | iso1+year,
  data = gspdat, cluster = "iso1"
))

summary(ptdspec7 <- feols(
  ipr_specific_prov ~ ptd_ma7 +
    log(gdppc_ppp+1) + log(gdp + 1) + gdp_growth_pcnt +
    democ + rta_signed | iso1+year,
  data = gspdat, cluster = "iso1"
))

summary(ptdsbst7 <- feols(
  ipr_scope_substantial_dummy ~ ptd_ma7 +
    log(gdppc_ppp+1) + log(gdp + 1) + gdp_growth_pcnt +
    democ + rta_signed | iso1+year,
  data = gspdat, cluster = "iso1"
))

```




## PTD analyses - results table
```{r}
etable(ptdtrips, ptdspec, ptdenf, ptdsbst,
       tex = T, depvar = T, digits = "r3", digits.stats = "r3")
```



```{r}
etable(ptdtrips5, ptdtrips7, ptdspec5, ptdspec7, ptdenf5, ptdenf7, ptdsbst5, ptdsbst7,
       tex = T, depvar = T, digits = "r3", digits.stats = "r3")
```




## PTD analyses - coefficient plot
```{r}
ptdmodels <- list(
  ptdtrips,
  ptdspec,
  ptdenf,
  ptdsbst
)

# Extract just the PTD coefficient + 95% CIs from each model
ptdcoef_df <- purrr::imap_dfr(ptdmodels, ~ {
  broom::tidy(.x, conf.int = TRUE) %>%        # get estimates + CIs
    filter(term == "ptd_ma") %>%                 # keep only PTD
    transmute(
      model      = .y,                        # model name
      estimate   = estimate,                  # point estimate
      conf.low   = conf.low,                  # lower 95% CI
      conf.high  = conf.high                  # upper 95% CI
    )
})

ptdcoef_df <- ptdcoef_df %>%
  mutate(
    model = factor(model,
      levels = c(1, 2, 3, 4),
      labels = c("TRIPS",
                 "Specific prov",
                 "Enforcement",
                 "Substantial"),
      ordered = T
    )
  )

# Make the coefficient plot
ggplot(ptdcoef_df, aes(x = estimate, y = model)) +
  geom_point() +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0, alpha=.5) +
  geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.5, color = "red", alpha=.5) +
  labs(
    x = "Coefficient estimates w/ 95% CIs",
    y = ""
  ) +
  xlim(c(-0.02,0.04)) +
  theme_minimal()
```



# Data cleansing
## US foreign aid
```{r}
usaid = read.csv("us_foreign_aid_country.csv")
usaid = usaid[!grepl("region", usaid$Country.Name, ignore.case = T),]
names(usaid)[names(usaid) == "Fiscal.Year"] = "year" # oct1-sep30 (year in which it ends)

desta = treaties_dyads
desta_detail = read.csv("desta_version_02_02.csv")
desta_ipr = desta_detail %>%
  select(number, ipr_specific_prov, ipr_pharma, ipr_trips_1994_dummy, ipr_scope_substantial_dummy, ipr_specific_enforcement_dummy)

desta = desta %>%
  left_join(desta_ipr, by = "number")

wto_dev = c("Albania", "Antigua", "Argentina", "Armenia",
            "Bahrain", "Barbados", "Belize", "Bolivia",
            "Botswana", "Brazil", "Brunei", "Verde",
            "Cameroon", "Chile", "China", "Taipei",
            "Colombia", "Congo", "Costa Rica", "Ivoire",
            "Cuba", "Dominica", "Dominican Republic", "Ecuador",
            "Egypt", "El Salvador", "Eswatini", "Fiji", "Gabon",
            "Georgia", "Ghana", "Grenada", "Guatemala", "Guyana",
            "Honduras", "Hong Kong", "India", "Indonesia",
            "Israel", "Jamaica", "Jordan", "Kazakhstan",
            "Kenya", "Korea", "Kuwait", "Kyrgyz", "Macao",
            "Malaysia", "Maldives", "Mauritius", "Mexico",
            "Moldova", "Mongolia", "Montenegro", "Morocco",
            "Namibia", "Nicaragua", "Nigeria", "Macedonia",
            "Oman", "Pakistan", "Panama", "Papua New Guinea",
            "Paraguay", "Peru", "Philippines", "Qatar",
            "Kitts", "Lucia", "Grenadines",
            "Samoa", "Saudi Arabia", "Seychelles", "Singapore",
            "South Africa", "Sri Lanka", "Suriname", "Tajikistan",
            "Thailand", "Tonga", "Trinidad", "Tunisia", "Turkey",
            "Ukraine", "United Arab Emirates", "Uruguay", "Vanuatu",
            "Venezuela", "Viet", "Zimbabwe") # source: tfadatabase.org

wto_ldc = c("Afghanistan", "Angola", "Bangladesh", "Benin",
            "Burkina Faso", "Burundi", "Cambodia", "Central African Republic",
            "Chad", "Congo", "Djibouti", "Gambia", "Guinea", "Bissau",
            "Haiti", "Laos", "Lesotho", "Liberia", "Madagascar",
            "Malawi", "Mali", "Mauritania", "Mozambique", "Myanmar",
            "Nepal", "Niger", "Rwanda", "Senegal", "Sierra Leone", "Solomon",
            "Tanzania", "Togo", "Uganda", "Yemen", "Zambia")

desta$devstat1 = ifelse(grepl(paste0(c(wto_dev, wto_ldc), collapse = "|"),
                              desta$country1, ignore.case = T), 1, 0)
desta$devstat2 = ifelse(grepl(paste0(c(wto_dev, wto_ldc), collapse = "|"),
                              desta$country2, ignore.case = T), 1, 0)

desta$devstat1 = replace(desta$devstat1, 
                         grepl("korea", desta$country1, ignore.case = T) 
                         & desta$year > 2019, 0)
desta$devstat2 = replace(desta$devstat2, 
                         grepl("korea", desta$country2, ignore.case = T) 
                         & desta$year > 2019, 0)
usdesta = desta[grepl("United States", desta$country1)|
                  grepl("United States", desta$country2),]
usdesta$iso1 = replace(usdesta$iso1, usdesta$iso1==840, 704)
usdesta$country1 = replace(usdesta$country1, usdesta$country1=="United States", "Vietnam")

usaid$devstat = ifelse(grepl(paste0(c(wto_dev, wto_ldc), collapse = "|"),
                              usaid$Country.Name, ignore.case = T), 1, 0)
usaid$devstat = replace(usaid$devstat, 
                         grepl("korea", usaid$Country.Name, ignore.case = T) 
                         & usaid$year > 2019, 0)
usaid = usaid[usaid$devstat==1,]
usaid$iso1 = countrycode(usaid$Country.Code, origin = "iso3c", destination = "iso3n")
usaid$year = as.numeric(usaid$year)
usaid = usaid[!is.na(usaid$year),]
usaid = usaid %>%
  left_join(usdesta, by = c("iso1", "year"))
usaid = usaid %>%
  select(-c(country1, country2, iso2))
usaid = usaid %>%
  select(-c(devstat, devstat1, devstat2))
usaid[is.na(usaid)] = 0


  
```


## OECD DAC
```{r}
dacdt = read.csv("OECD_DAC.csv")
dacdt = dacdt %>%
  mutate(year = TIME_PERIOD,
         iso1 = countrycode(RECIPIENT, "iso3c", "iso3n"),
         donor = DONOR,
         ODAusd = OBS_VALUE*10^6) %>%
  filter(!is.na(iso1)) %>%
  select(iso1, donor, year, ODAusd)

us_dac = dacdt %>%
  filter(donor == "USA") %>%
  select(-donor)

eu_dac = dacdt %>%
  filter(donor == "DACEU") %>%
  select(-donor)

us_dac = us_dac %>%
  left_join(us_desta, by = c("iso1", "year"))

eu_dac = eu_dac %>%
  left_join(eu_desta, by = c("iso1", "year")) %>%
  dplyr::group_by(iso1, year) %>%
  mutate(ODAusd = sum(ODAusd)) %>%
  distinct(iso1, year, .keep_all = T)

```




## covariates - trade stats
```{r}
## import and export
tradeflow_us = read.csv("dots_us.csv")
tradeflow_us = tradeflow_us %>% select(
  Counterpart.Country.Name,
  Counterpart.Country.Code,
  Time.Period,
  Goods..Value.of.Imports..Cost..Insurance..Freight..CIF...US.Dollars..TMG_CIF_USD.,
  Goods..Value.of.Exports..Free.on.board..FOB...US.Dollars..TXG_FOB_USD.
)
colnames(tradeflow_us) = c("partner", "cc_imf", "year", "import", "export")
tradeflow_us = tradeflow_us %>% mutate(
  iso1 = countrycode(cc_imf, "imf", "iso3n")
)

tradeflow_us = tradeflow_us[!is.na(tradeflow_us$iso1),]
tradeflow_us = tradeflow_us %>% select(
  -c(partner, cc_imf)
)

tradeflow_eu = read.csv("dots_eu.csv")
tradeflow_eu = tradeflow_eu %>% select(
  Counterpart.Country.Name,
  Counterpart.Country.Code,
  Time.Period,
  Goods..Value.of.Imports..Cost..Insurance..Freight..CIF...US.Dollars..TMG_CIF_USD.,
  Goods..Value.of.Exports..Free.on.board..FOB...US.Dollars..TXG_FOB_USD.
)
colnames(tradeflow_eu) = c("partner", "cc_imf", "year", "import", "export")
tradeflow_eu = tradeflow_eu %>% mutate(
  iso1 = countrycode(cc_imf, "imf", "iso3n")
)

tradeflow_eu = tradeflow_eu[!is.na(tradeflow_eu$iso1),]
tradeflow_eu = tradeflow_eu %>% select(
  -c(partner, cc_imf)
)

usaid = usaid %>%
  left_join(
    tradeflow_us,
    by = c("iso1", "year")
  )

us_dac = us_dac %>%
  left_join(
    tradeflow_us,
    by = c("iso1", "year")
  )

eu_dac = eu_dac %>%
  left_join(
    tradeflow_eu,
    by = c("iso1", "year")
  )

## GDP pc
gdppc = WDI(
  country = "all",
  indicator = "NY.GDP.PCAP.PP.KD",
  start = 1960,
  end = 2023,
  extra = T
)

gdppc = gdppc %>%
  mutate(
    iso1 = countrycode(iso3c, "iso3c", "iso3n") 
  )

gdppc = gdppc[!is.na(gdppc$iso1),]

gdppc = gdppc %>% select(
  iso1,
  year,
  NY.GDP.PCAP.PP.KD
)

names(gdppc)[names(gdppc) == "NY.GDP.PCAP.PP.KD"] = "gdppc_ppp"

usaid = usaid %>%
  left_join(
    gdppc,
    by = c("iso1", "year")
  )

us_dac = us_dac %>%
  left_join(
    gdppc,
    by = c("iso1", "year")
  )

eu_dac = eu_dac %>%
  left_join(
    gdppc,
    by = c("iso1", "year")
  )

## GDP constant / GDP growth
gdpcons = WDI(
  country = "all",
  indicator = "NY.GDP.MKTP.PP.KD", #GDP PPP, constant 2021 USD
  start = 1960,
  end = 2023,
  extra = T
)

gdpgrowth = WDI(
  country = "all",
  indicator = "NY.GDP.MKTP.KD.ZG",
  start = 1960,
  end = 2023,
  extra = T
)

gdpcons = gdpcons %>% select(iso3c, year, NY.GDP.MKTP.PP.KD)
gdpgrowth = gdpgrowth %>% select(iso3c, year, NY.GDP.MKTP.KD.ZG)

gdpcons$iso1 = countrycode(gdpcons$iso3c, "iso3c", "iso3n")
gdpgrowth$iso1 = countrycode(gdpgrowth$iso3c, "iso3c", "iso3n")

usaid = usaid %>%
  left_join(
    gdpcons,
    by = c("iso1", "year")
  )

us_dac = us_dac %>%
  left_join(
    gdpcons,
    by = c("iso1", "year")
  )

eu_dac = eu_dac %>%
  left_join(
    gdpcons,
    by = c("iso1", "year")
  )

usaid = usaid %>%
  left_join(
    gdpgrowth,
    by = c("iso1", "year")
  )

us_dac = us_dac %>%
  left_join(
    gdpgrowth,
    by = c("iso1", "year")
  )

eu_dac = eu_dac %>%
  left_join(
    gdpgrowth,
    by = c("iso1", "year")
  )

usaid = usaid %>% select(-iso3c.x, -iso3c.y)
names(usaid)[names(usaid) == "NY.GDP.MKTP.PP.KD"] = "gdp"
names(usaid)[names(usaid) == "NY.GDP.MKTP.KD.ZG"] = "gdp_growth_pcnt"

names(us_dac)[names(us_dac) == "NY.GDP.MKTP.PP.KD"] = "gdp"
names(us_dac)[names(us_dac) == "NY.GDP.MKTP.KD.ZG"] = "gdp_growth_pcnt"

names(eu_dac)[names(eu_dac) == "NY.GDP.MKTP.PP.KD"] = "gdp"
names(eu_dac)[names(eu_dac) == "NY.GDP.MKTP.KD.ZG"] = "gdp_growth_pcnt"

## RTAs signed
desta = desta %>% mutate(
  which_cty = paste(country1, country2, sep = " ")
)

desta_dev = desta %>% filter(
  grepl(paste0(c(wto_dev, wto_ldc), collapse = "|"),
        which_cty))

desta_dev = desta_dev %>% mutate(
    temp = ifelse(devstat1 == 0, country1, country2),
    temp2 = ifelse(devstat1 == 0, iso1, iso2),
    iso1 = ifelse(devstat1 == 0, iso2, iso1),
    iso2 = ifelse(devstat1== 0, temp2, iso2),
    country1 = ifelse(devstat1 == 0, country2, country1),
    country2 = ifelse(devstat1 == 0, temp, country2)
  ) %>%
  select(-c(temp, temp2))

desta_dev = desta_dev %>% distinct(country1, name, .keep_all = T)

desta_dev = desta_dev %>% filter(entry_type == "base_treaty")

desta_dev = desta_dev %>%
  group_by(country1, year) %>%
  mutate(rta_signed = n())
desta_dev = desta_dev %>% distinct(iso1, year, .keep_all = T)
desta_dev = desta_dev %>% ungroup() %>% select(iso1, year, rta_signed)

usaid = usaid %>%
  left_join(
    desta_dev,
    by = c("iso1", "year")
  )

us_dac = us_dac %>%
  left_join(
    desta_dev,
    by = c("iso1", "year")
  )

eu_dac = eu_dac %>%
  left_join(
    desta_dev,
    by = c("iso1", "year")
  )

usaid$rta_signed = replace(usaid$rta_signed, is.na(usaid$rta_signed), 0)
us_dac = us_dac %>% replace_na(list(rta_signed = 0))
eu_dac = eu_dac %>% replace_na(list(rta_signed = 0))

## Polity
polity = read_spss("p5v2018.sav")
polity = polity %>% filter(year >= 1990)
polity = polity %>% mutate(
  iso1 = countrycode(ccode, "cown", "iso3n")
)
polity = polity %>% select(iso1, year, polity2)

usaid = usaid %>%
  left_join(
    polity,
    by = c("iso1", "year")
  )

us_dac = us_dac %>%
  left_join(
    polity,
    by = c("iso1", "year")
  )

eu_dac = eu_dac %>%
  left_join(
    polity,
    by = c("iso1", "year")
  )

usaid = usaid %>%
  group_by(Country.Name) %>%
  fill(polity2, .direction = "updown")
usaid = usaid %>% mutate(
  democ = ifelse(polity2 > 0, 1, 0)
)

us_dac = us_dac %>%
  group_by(iso1) %>%
  fill(polity2, .direction = "updown") %>%
  ungroup() %>%
  mutate(
  democ = ifelse(polity2 > 0, 1, 0)
)

eu_dac = eu_dac %>%
  group_by(iso1) %>%
  fill(polity2, .direction = "updown") %>% 
  ungroup() %>%
  mutate(
  democ = ifelse(polity2 > 0, 1, 0)
)

usaid = usaid %>%
  mutate(
    ldc = ifelse(Income.Group.Acronym == "LIC", 1, 0)
  )

usaid_disb = usaid[usaid$Transaction.Type.Name == "Disbursements",]
usaid_obl = usaid[usaid$Transaction.Type.Name == "Obligations",]
usaid_app = usaid[usaid$Transaction.Type.Name == "Appropriated and Planned",]


```




# US aid prelim analysis
## TRIPS
```{r}
# constant_amount is the US aid in current USD
# ldc multicollinearity, removed
summary(feols(
  ipr_trips_1994_dummy ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + ldc|iso1+year, 
  data = usaid_app, cluster = "iso1"))
summary(feols(
  ipr_trips_1994_dummy ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + ldc|iso1+year, 
  data = usaid_obl, cluster = "iso1"))
summary(feols(
  ipr_trips_1994_dummy ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + ldc|iso1+year, 
  data = usaid_disb, cluster = "iso1"))
```


## Specific provisions
```{r}
summary(feols(
  ipr_specific_enforcement_dummy ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1)|iso1+year, 
  data = usaid_app, cluster = "iso1"))
summary(feols(
  ipr_specific_enforcement_dummy ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1)|iso1+year, 
  data = usaid_obl, cluster = "iso1"))
summary(feols(
  ipr_specific_enforcement_dummy ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1)|iso1+year, 
  data = usaid_disb, cluster = "iso1"))
```



## Pharmaceuticals
```{r}
summary(feols(
  ipr_pharma ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1)|iso1+year, 
  data = usaid_app, cluster = "iso1"))
summary(feols(
  ipr_pharma ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1)|iso1+year, 
  data = usaid_obl, cluster = "iso1"))
summary(feols(
  ipr_pharma ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1)|iso1+year, 
  data = usaid_disb, cluster = "iso1"))
```


# US aid analyses w/ covariates
## TRIPS

```{r}
summary(feols(
  ipr_trips_1994_dummy ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ + ldc|iso1+year, 
  data = usaid_app, cluster = "iso1"))
summary(feols(
  ipr_trips_1994_dummy ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ + ldc|iso1+year, 
  data = usaid_obl, cluster = "iso1"))
summary(aidtrips <- feols(
  ipr_trips_1994_dummy ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ + ldc|iso1+year, 
  data = usaid_disb, cluster = "iso1"))
```



```{r}
eutemp = eu_dac %>%
  dplyr::group_by(iso1) %>%
  dplyr::arrange(iso1, year) %>%
  mutate(logODA = sign(ODAusd)*log1p(abs(ODAusd))) %>%
  filter(!is.na(logODA)) %>%
  mutate(llogODA = lag(logODA)) %>%
  ungroup() %>%
  mutate(across(starts_with("ipr_"), ~replace_na(.x, 0)))

summary(eutrips <- feols(
  ipr_trips_1994_dummy ~ llogODA + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year, 
  data = eutemp, cluster = "iso1"))

```



## Specific provisions

```{r}
summary(feols(
  ipr_specific_enforcement_dummy ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ + ldc|iso1+year, 
  data = usaid_app, cluster = "iso1"))
summary(feols(
  ipr_specific_enforcement_dummy ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ + ldc|iso1+year, 
  data = usaid_obl, cluster = "iso1"))
summary(aidenf <- feols(
  ipr_specific_enforcement_dummy ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ + ldc|iso1+year, 
  data = usaid_disb, cluster = "iso1"))
```



```{r}
summary(euenf <- feols(
  ipr_specific_enforcement_dummy ~ llogODA + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year, 
  data = eutemp, cluster = "iso1"))

```



```{r}
summary(feols(
  ipr_specific_prov ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ + ldc|iso1+year, 
  data = usaid_app, cluster = "iso1"))
summary(feols(
  ipr_specific_prov ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ + ldc|iso1+year, 
  data = usaid_obl, cluster = "iso1"))
summary(aidspec <- feols(
  ipr_specific_prov ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ + ldc|iso1+year, 
  data = usaid_disb, cluster = "iso1"))
```



```{r}
summary(euspec <- feols(
  ipr_specific_prov ~ llogODA + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year, 
  data = eutemp, cluster = "iso1"))
```




```{r}
summary(feols(
  ipr_scope_substantial_dummy ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ + ldc|iso1+year, 
  data = usaid_app, cluster = "iso1"))
summary(feols(
  ipr_scope_substantial_dummy ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ + ldc|iso1+year, 
  data = usaid_obl, cluster = "iso1"))
summary(aidsbst <- feols(
  ipr_scope_substantial_dummy ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ + ldc|iso1+year, 
  data = usaid_disb, cluster = "iso1"))
```


```{r}
summary(eusbst <- feols(
  ipr_scope_substantial_dummy ~ llogODA + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year, 
  data = eutemp, cluster = "iso1"))
```




## Pharmaceuticals

```{r}
summary(feols(
  ipr_pharma ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ + ldc|iso1+year, 
  data = usaid_app, cluster = "iso1"))
summary(feols(
  ipr_pharma ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ + ldc|iso1+year, 
  data = usaid_obl, cluster = "iso1"))
summary(feols(
  ipr_pharma ~ log(constant_amount+1) + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ + ldc|iso1+year, 
  data = usaid_disb, cluster = "iso1"))
```



```{r}
summary(feols(
  ipr_pharma ~ llogODA + 
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year, 
  data = eutemp, cluster = "iso1"))
```



## US aid analyses - results table
```{r}
etable(aidtrips, aidspec, aidenf, aidsbst,
       tex = T, depvar = T, digits = "r3", digits.stats = "r3")
```



## US aid analyses - coefficient plot
```{r}
usaidmodels <- list(
  aidtrips,
  aidspec,
  aidenf,
  aidsbst
)

# Extract just the PTD coefficient + 95% CIs from each model
usaidcoef_df <- purrr::imap_dfr(usaidmodels, ~ {
  broom::tidy(.x, conf.int = TRUE) %>%        # get estimates + CIs
    filter(term == "log(constant_amount + 1)") %>%                 # keep only aid
    transmute(
      model      = .y,                        # model name
      estimate   = estimate,                  # point estimate
      conf.low   = conf.low,                  # lower 95% CI
      conf.high  = conf.high                  # upper 95% CI
    )
})

usaidcoef_df <- usaidcoef_df %>%
  mutate(
    model = factor(model,
      levels = c(1, 2, 3, 4),
      labels = c("TRIPS",
                 "Specific prov",
                 "Enforcement",
                 "Substantial"),
      ordered = T
    )
  )

# Make the coefficient plot
ggplot(usaidcoef_df, aes(x = estimate, y = model)) +
  geom_point() +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0, alpha=.5) +
  geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.5, color = "red", alpha =.5) +
  labs(
    x = "Coefficient estimates w/ 95% CIs",
    y = ""
  ) +
  xlim(c(-0.003, 0.00875)) + 
  theme_minimal()
```



## EU aid analyses - results table
```{r}
etable(eutrips, euspec, euenf, eusbst,
       tex = T, depvar = T, digits = "r3", digits.stats = "r3")
```



## EU aid analyses - coefficient plot
```{r}
euaidmodels <- list(
  eutrips,
  euspec,
  euenf,
  eusbst
)

# Extract just the PTD coefficient + 95% CIs from each model
euaidcoef_df <- purrr::imap_dfr(euaidmodels, ~ {
  broom::tidy(.x, conf.int = TRUE) %>%        # get estimates + CIs
    filter(term == "llogODA") %>%                 # keep only aid
    transmute(
      model      = .y,                        # model name
      estimate   = estimate,                  # point estimate
      conf.low   = conf.low,                  # lower 95% CI
      conf.high  = conf.high                  # upper 95% CI
    )
})

euaidcoef_df <- euaidcoef_df %>%
  mutate(
    model = factor(model,
      levels = c(1, 2, 3, 4),
      labels = c("TRIPS",
                 "Specific prov",
                 "Enforcement",
                 "Substantial"),
      ordered = T
    )
  )

# Make the coefficient plot
ggplot(euaidcoef_df, aes(x = estimate, y = model)) +
  geom_point() +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0, alpha=.5) +
  geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.5, color = "red", alpha=.5) +
  labs(
    x = "Coefficient estimates w/ 95% CIs",
    y = ""
  ) +
  theme_minimal()
```



# Summary stats
```{r, warning=FALSE}
summary(gspdat%>%
          ungroup() %>%
          select(ptd_ma,
                 gdppc_ppp, gdp, gdp_growth_pcnt,
                 rta_signed, democ,
                 ipr_trips_1994_dummy,
                 ipr_specific_prov,
                 ipr_specific_enforcement_dummy,
                 ipr_scope_substantial_dummy)%>%
          mutate(lgdppc=log(gdppc_ppp+1),
                 lgdp=log(gdp+1)) %>%
          select(ptd_ma, lgdppc, lgdp,
                 gdp_growth_pcnt, rta_signed, democ,
                 ipr_trips_1994_dummy,
                 ipr_specific_prov,
                 ipr_specific_enforcement_dummy,
                 ipr_scope_substantial_dummy)
          )

gspdat%>%
  ungroup() %>%
  select(ptd_ma, gdppc_ppp, gdp, 
         gdp_growth_pcnt,rta_signed, democ,
         ipr_trips_1994_dummy,
                 ipr_specific_prov,
                 ipr_specific_enforcement_dummy,
                 ipr_scope_substantial_dummy) %>% 
  mutate(lgdppc=log(gdppc_ppp+1),
                 lgdp=log(gdp+1)) %>% 
  select(ptd_ma, lgdppc, lgdp,
                 gdp_growth_pcnt, rta_signed, democ,
         ipr_trips_1994_dummy,
                 ipr_specific_prov,
                 ipr_specific_enforcement_dummy,
                 ipr_scope_substantial_dummy) %>% 
  dplyr::summarise(across(c(ptd_ma, lgdppc, lgdp, 
                            gdp_growth_pcnt, rta_signed, democ,
                            ipr_trips_1994_dummy,
                 ipr_specific_prov,
                 ipr_specific_enforcement_dummy,
                 ipr_scope_substantial_dummy),.fns = sd,na.rm = T))

summary(usaid_disb%>%
          ungroup() %>%
          select(constant_amount, import, export)%>%
          mutate(logaid=log(constant_amount+1),
                 limport=log(import+1),
                 lexport=log(export+1)) %>%
          select(logaid, limport, lexport))

usaid_disb%>%
          ungroup() %>%
          select(constant_amount, import, export)%>%
          mutate(logaid=log(constant_amount+1),
                 limport=log(import+1),
                 lexport=log(export+1)) %>%
          select(logaid, limport, lexport) %>%
  dplyr::summarise(across(c(logaid, limport, lexport), .fns = sd, na.rm=T))

summary(eutemp %>%
          ungroup() %>%
          select(llogODA, import, export) %>%
          mutate(limport=log(import+1),
                 lexport=log(export+1)) %>%
          select(llogODA, limport, lexport))

eutemp %>%
  ungroup() %>%
  select(llogODA, import, export) %>%
  mutate(limport=log(import+1),
         lexport=log(export+1)) %>%
  select(llogODA, limport, lexport) %>%
  dplyr::summarise(across(c(llogODA, limport, lexport), .fns = sd, na.rm=T))

```



# Endogenous aid - 2SLS
## Mortality data prep
```{r}
mort = read.csv("infec_mortality.csv")
mort = mort %>%
  mutate(deaths=log1p(val),
         iso3c = countryname(location_name, destination = "iso3c")) 
mort = mort %>%
  group_by(iso3c) %>%
  dplyr::arrange(iso3c, year) %>%
  mutate(
    mdeath = mean(val, na.rm=T),
    sdeath = sd(val, na.rm = T),
    std_deaths = (val-mdeath)/sdeath)

mort$iso3c = replace(mort$iso3c, is.na(mort$iso3c)&grepl("lebanese", mort$location_name, ignore.case = T), "LBN")
mort$iso3c = replace(mort$iso3c, is.na(mort$iso3c)&grepl("portugese", mort$location_name, ignore.case = T), "PRT")
mort$iso3c = replace(mort$iso3c, is.na(mort$iso3c)&grepl("guyana", mort$location_name, ignore.case = T), "GUY")
mort = mort %>%
  dplyr::mutate(iso1 = countrycode(iso3c, "iso3c", "iso3n")) %>%
  dplyr::ungroup() %>%
  select(iso1, year, deaths, std_deaths)

disaster_mort = read.csv("disaster_mortality.csv")
disaster_mort = disaster_mort %>%
  mutate(dist_deaths=log1p(val),
         iso3c = countryname(location_name, destination = "iso3c")) 
disaster_mort = disaster_mort %>%
  group_by(iso3c) %>%
  dplyr::arrange(iso3c, year) %>%
  mutate(
    mdeath = mean(val, na.rm=T),
    sdeath = sd(val, na.rm = T),
    std_dist_deaths = (val-mdeath)/sdeath)
disaster_mort$iso3c = replace(disaster_mort$iso3c, is.na(disaster_mort$iso3c)&grepl("lebanese", disaster_mort$location_name, ignore.case = T), "LBN")
disaster_mort$iso3c = replace(disaster_mort$iso3c, is.na(disaster_mort$iso3c)&grepl("portugese", disaster_mort$location_name, ignore.case = T), "PRT")
disaster_mort$iso3c = replace(disaster_mort$iso3c, is.na(disaster_mort$iso3c)&grepl("guyana", disaster_mort$location_name, ignore.case = T), "GUY")
disaster_mort = disaster_mort %>%
  dplyr::mutate(iso1 = countrycode(iso3c, "iso3c", "iso3n")) %>%
  dplyr::ungroup() %>%
  select(iso1, year, dist_deaths, std_dist_deaths)

unga_align = read.csv("UNGA_alignment.csv")
unga_align = unga_align %>%
  filter(year>=1990) %>%
  mutate(iso1 = countrycode(member_state_id, origin = "cown", destination = "iso3n"),
         unga_impvote = S_USA_important) %>%
  select(iso1, year, unga_impvote)

usaid_disb = usaid_disb %>%
  left_join(unga_align, by=c("iso1", "year"))


usaid_disb2 = usaid_disb %>%
  left_join(mort, by=c("iso1", "year")) %>%
  left_join(disaster_mort, by=c("iso1", "year")) %>%
  replace_na(list(deaths=0, std_deaths=0, dist_deaths=0, std_dist_deaths=0)) %>%
  dplyr::mutate(logaid = sign(constant_amount)*log1p(abs(constant_amount))) %>%
  dplyr::group_by(iso1) %>%
  dplyr::arrange(iso1, year) %>%
  dplyr::mutate(ldeaths = lag(deaths),
                lstd_deaths = lag(std_deaths),
                ldist_deaths = lag(dist_deaths),
                llogaid = lag(log1p(constant_amount)),
                limpvote = lag(unga_impvote))


eutemp2 = eutemp %>%
  left_join(mort, by=c("iso1", "year")) %>%
  replace_na(list(deaths=0, std_deaths=0)) %>%
  group_by(iso1) %>%
  dplyr::mutate(ldeaths = lag(deaths,2),
                lstd_deaths = lag(std_deaths,2))

```



## UNSC data prep
```{r}
unsc_membership <- list(
  "Algeria" = 2004:2005,
  "Angola" = c(2003:2004, 2015:2016),
  "Argentina" = c(1995, 1999:2000, 2005:2006, 2013:2014),
  "Bahrain" = 1998:1999,
  "Bangladesh" = 2000:2001,
  "Bolivia" = 2017:2018,
  "Botswana" = 1995:1996,
  "Brazil" = c(1998:1999, 2004:2005, 2010:2011, 2022),
  "Bulgaria" = 2002:2003,
  "Cameroon" = 2002:2003,
  "Chile" = c(1996:1997, 2003:2004, 2014:2015),
  "Colombia" = c(2001:2002, 2011:2012),
  "Costa Rica" = c(1997:1998, 2008:2009),
  "Cote D’ivoire" = 2018:2019,
  "Ecuador" = 2017:2018,
  "Egypt" = c(1996:1997, 2016:2017),
  "Ethiopia" = 2017:2018,
  "Ghana" = 2006:2007,
  "Guatemala" = 2012:2013,
  "India" = c(2011:2012),
  "Indonesia" = c(1995:1996, 2007:2008),
  "Jordan" = 2014:2015,
  "Kenya" = 1997:1998,
  "Korea" = c(1996:1997, 2013:2014),  
  "Malaysia" = 2015:2016,
  "Mexico" = c(2002:2003, 2009:2010),
  "Morocco" = c(2012:2013),
  "Namibia" = 1999:2000,
  "New Zealand" = c(2015:2016),
  "Nigeria" = c(1995, 2010:2011),
  "Pakistan" = c(2012:2013),
  "Panama" = 2007:2008,
  "Peru" = c(2006:2007, 2018:2019),
  "Philippines" = 2004:2005,
  "Poland" = c(1996:1997, 2018:2019),
  "Romania" = c(2004:2005),
  "South Africa" = c(2007:2008, 2011:2012),
  "Tunisia" = 2000:2001,
  "Turkey" = 2009:2010,
  "Ukraine" = c(2000:2001, 2016:2017),
  "United Arab Emirates" = 2022,
  "Uruguay" = 2016:2017,
  "Venezuela" = 2015:2016,
  "Vietnam" = 2008:2009
)

unsc_df <- tibble(
  country = rep(names(unsc_membership), lengths(unsc_membership)),
  year = unlist(unsc_membership)
) %>%
  mutate(unsc = 1,
         iso1 = countrycode(country, "country.name", "iso3n"))

usaid_disb2 = usaid_obl %>%
  left_join(unsc_df, by=c("iso1", "year")) %>%
  replace_na(list(unsc=0)) %>%
  dplyr::mutate(logaid = sign(constant_amount)*log1p(abs(constant_amount))) %>%
  dplyr::group_by(iso1) %>%
  dplyr::arrange(iso1, year) %>%
  dplyr::mutate(llogaid = lag(logaid),
                lunsc = lag(unsc))

```



## US aid- mortality only
```{r}
summary(aidtrips_2sls0 <- feols(
  ipr_trips_1994_dummy ~
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year|
    log1p(constant_amount) ~ ldeaths, 
  data = usaid_disb2, cluster = "iso1"))

summary(aidspec_2sls0 <- feols(
  ipr_specific_prov ~
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year|
    log1p(constant_amount) ~ ldeaths, 
  data = usaid_disb2, cluster = "iso1"))

summary(aidenf_2sls0 <- feols(
  ipr_specific_enforcement_dummy ~
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year|
    log1p(constant_amount) ~ ldeaths, 
  data = usaid_disb2, cluster = "iso1"))

summary(aidsbst_2sls0 <- feols(
  ipr_scope_substantial_dummy ~
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year|
    log1p(constant_amount) ~ ldeaths, 
  data = usaid_disb2, cluster = "iso1"))
```


## US aid - mortality + imp votes
```{r}
summary(aidtrips_2sls <- feols(
  ipr_trips_1994_dummy ~
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year|
    log1p(constant_amount) ~ ldeaths + limpvote, 
  data = usaid_disb2, cluster = "iso1"))

summary(aidspec_2sls <- feols(
  ipr_specific_prov ~
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year|
    log1p(constant_amount) ~ ldeaths + limpvote, 
  data = usaid_disb2, cluster = "iso1"))

summary(aidenf_2sls <- feols(
  ipr_specific_enforcement_dummy ~
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year|
    log1p(constant_amount) ~ ldeaths + limpvote, 
  data = usaid_disb2, cluster = "iso1"))

summary(aidsbst_2sls <- feols(
  ipr_scope_substantial_dummy ~
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year|
    log1p(constant_amount) ~ ldeaths + limpvote, 
  data = usaid_disb2, cluster = "iso1"))
```



```{r}
etable(aidtrips_2sls0, aidtrips_2sls, aidspec_2sls0, aidspec_2sls,
       aidenf_2sls0, aidenf_2sls, aidsbst_2sls0, aidsbst_2sls,
       tex = T, depvar = T, digits = "r3", digits.stats = "r3")
```



## coefplot
```{r}
usaid2slsmodels <- list(
  aidtrips_2sls0, aidtrips_2sls, aidspec_2sls0, aidspec_2sls,
       aidenf_2sls0, aidenf_2sls, aidsbst_2sls0, aidsbst_2sls
)

# Extract just the PTD coefficient + 95% CIs from each model
usaid2slscoef_df <- purrr::imap_dfr(usaid2slsmodels, ~ {
  broom::tidy(.x, conf.int = TRUE) %>%        # get estimates + CIs
    filter(term == "fit_log1p(constant_amount)") %>%                 # keep only aid
    transmute(
      model      = .y,                        # model name
      estimate   = estimate,                  # point estimate
      conf.low   = conf.low,                  # lower 95% CI
      conf.high  = conf.high                  # upper 95% CI
    )
})

usaid2slscoef_df <- usaid2slscoef_df %>%
  mutate(
    model = factor(model,
      levels = c(1, 2, 3, 4, 5, 6, 7, 8),
      labels = c("TRIPS", "TRIPS",
                 "Specific prov", "Specific prov",
                 "Enforcement", "Enforcement",
                 "Substantial", "Substantial"),
      ordered = T
    )
  )
usaid2slscoef_df$IV = rep(c("Mortality", "Mortality + Imp. votes"), 4)

# Make the coefficient plot
ggplot(usaid2slscoef_df, aes(x = estimate, y = model, color = IV)) +
  geom_point(position=position_dodge(width=0.5)) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), position=position_dodge(width=0.5), width = 0, alpha=.5) +
  scale_color_manual(values=c("grey55", "black"), name = "IV") +
  geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.5, color = "red", alpha =.5) +
  labs(
    x = "Coefficient estimates w/ 95% CIs",
    y = ""
  ) +
  theme_minimal() + 
  xlim(c(-0.05, 0.1))
```



## EU aid
```{r}
summary(eutrips_2sls <- feols(
  ipr_trips_1994_dummy ~
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year|
    llogODA ~ ldeaths, 
  data = eutemp2, cluster = "iso1"))

summary(euspec_2sls <- feols(
  ipr_specific_prov ~
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year|
    llogODA ~ ldeaths, 
  data = eutemp2, cluster = "iso1"))

summary(euenf_2sls <- feols(
  ipr_specific_enforcement_dummy ~
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year|
    llogODA ~ ldeaths, 
  data = eutemp2, cluster = "iso1"))

summary(eusbst_2sls <- feols(
  ipr_scope_substantial_dummy ~
    log(import+1) + log(export+1) + log(gdppc_ppp+1) + 
    log(gdp+1) + gdp_growth_pcnt + rta_signed + democ|iso1+year|
    llogODA ~ ldeaths, 
  data = eutemp2, cluster = "iso1"))
```


```{r}
etable(eutrips_2sls, euspec_2sls,
       euenf_2sls, eusbst_2sls,
       tex = T, depvar = T, digits = "r3", digits.stats = "r3")
```



## coefplot
```{r}
euaid2slsmodels <- list(
  eutrips_2sls,
  euspec_2sls,
  euenf_2sls,
  eusbst_2sls
)

# Extract just the PTD coefficient + 95% CIs from each model
euaid2slscoef_df <- purrr::imap_dfr(euaid2slsmodels, ~ {
  broom::tidy(.x, conf.int = TRUE) %>%        # get estimates + CIs
    filter(term == "fit_llogODA") %>%                 # keep only aid
    transmute(
      model      = .y,                        # model name
      estimate   = estimate,                  # point estimate
      conf.low   = conf.low,                  # lower 95% CI
      conf.high  = conf.high                  # upper 95% CI
    )
})

euaid2slscoef_df <- euaid2slscoef_df %>%
  mutate(
    model = factor(model,
      levels = c(1, 2, 3, 4),
      labels = c("TRIPS",
                 "Specific prov",
                 "Enforcement",
                 "Substantial"),
      ordered = T
    )
  )

# Make the coefficient plot
ggplot(euaid2slscoef_df, aes(x = estimate, y = model)) +
  geom_point() +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0, alpha=.5) +
  geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.5, color = "red", alpha=.5) +
  labs(
    x = "Coefficient estimates w/ 95% CIs",
    y = ""
  ) +
  theme_minimal() +
  xlim(c(-0.4, 0.4))
```





